library(dplyr)
library(ggplot2)
library(vtable)
library(fixest)
library(margins)
library(clubSandwich)
library(multiwayvcov)
library(lmtest)
library(glm.predict)

setwd("Directory")

load(file='Rule.RCL.Ideology.rda')
load(file='Rule.CL.Ideology.rda')


# Figure 1
# Note: As 'margins' and 'glm.predict' packages do not support 'fixest' objects,
# we calculate average marginal effects and predicted probability estimates based on the results from glm. 
# The differences in standard errors between feglm and glm are negligible and 
# do not change the significance of the variables. 
Baseline.glm <- glm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                      Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                      Internal.Joint + Emergency + Contract + 
                      Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                      Unalign.President + FirstYear + LastYear +  
                      Congressional.Opposion + Statutory.Deadline + 
                      CourtCaseAPA.log + Judicial.Deadline +
                      Ideo.Mismatch.norm.RCL + 
                      Agency_Bureau + President, 
                    data=Rule.RCL.Ideology, binomial(link="logit"))
Baseline.glm.cl <- cluster.vcov(Baseline.glm, ~Agency)
Baseline.glm.coef <- broom::tidy(coeftest(Baseline.glm, Baseline.glm.cl))

Baseline.AME <- margins(Baseline.glm, 
                        vcov = cluster.vcov(Baseline.glm, ~Agency),
                        variables = c("Rule.Novelty","Routine","Word.Count.log",
                                      "Stringency4.log","Econ.Sig","SL.Gov",
                                      "External.Joint","Internal.Joint","Emergency","Contract",
                                      "Tot.Emp.log","HQ.Pro.per","Rule.Burden.log",
                                      "Unalign.President","FirstYear","LastYear",
                                      "Congressional.Opposion","Statutory.Deadline",
                                      "CourtCaseAPA.log","Judicial.Deadline","Ideo.Mismatch.norm.RCL"),
                        type = "response") %>%
  summary()

Baseline.AME.changeSD <- margins(Baseline.glm, 
                                 vcov = cluster.vcov(Baseline.glm, ~Agency),
                                 variables = c("Word.Count.log",
                                               "Stringency4.log",
                                               "Tot.Emp.log","HQ.Pro.per","Rule.Burden.log",
                                               "Congressional.Opposion",
                                               "CourtCaseAPA.log","Ideo.Mismatch.norm.RCL"),
                                 change = "sd",
                                 type = "response") %>%
  summary()

Baseline.AME.all <- Baseline.AME %>%
  subset(factor!="Word.Count.log" & factor!="Stringency4.log" & factor!="Tot.Emp.log" &
         factor!="HQ.Pro.per" & factor!="Rule.Burden.log" & factor!="Congressional.Opposion" &
         factor!="CourtCaseAPA.log" & factor!="Ideo.Mismatch.norm.RCL") %>%
  rbind(Baseline.AME.changeSD) %>%
  mutate(Var=ifelse(factor=="Rule.Novelty1","New Rule",
             ifelse(factor=="Routine1","Routine and Frequent Rule",
             ifelse(factor=="Word.Count.log","*Rule Complexity",
             ifelse(factor=="Stringency4.log","*Rule Stringency",      
             ifelse(factor=="Econ.Sig1","Economic Significance",
             ifelse(factor=="SL.Gov1","Other Levels of Government",  
             ifelse(factor=="External.Joint1","Joint Rulemaking - External",  
             ifelse(factor=="Internal.Joint1","Joint Rulemaking - Internal",  
             ifelse(factor=="Emergency1","Emergency Rule",
             ifelse(factor=="Contract1","Contract Rule",
             ifelse(factor=="HQ.Pro.per","*Professionalism",
             ifelse(factor=="Tot.Emp.log","*Size (Employee)",
             ifelse(factor=="Rule.Burden.log","*Rulemaking Burden",
             ifelse(factor=="Unalign.President1","Unaligned President",
             ifelse(factor=="FirstYear1","Presidents’ First Year",
             ifelse(factor=="LastYear1","Presidents’ Last Year",
             ifelse(factor=="Congressional.Opposion","*Congressional Opposition",
             ifelse(factor=="Statutory.Deadline1","Statutory Deadline",
             ifelse(factor=="CourtCaseAPA.log","*Court Cases (APA)",
             ifelse(factor=="Judicial.Deadline1","Judicial Deadline",
             ifelse(factor=="Ideo.Mismatch.norm.RCL","*Interest Group Opposition","XXX")))))))))))))))))))))) %>%
  mutate(Var=factor(Var,
                    levels=c("New Rule","Routine and Frequent Rule","*Rule Complexity","*Rule Stringency",
                             "Economic Significance","Other Levels of Government","Joint Rulemaking - External",
                             "Joint Rulemaking - Internal","Emergency Rule","Contract Rule",
                             "*Professionalism","*Size (Employee)","*Rulemaking Burden",
                             "Unaligned President","Presidents’ First Year","Presidents’ Last Year",
                             "*Congressional Opposition","Statutory Deadline",
                             "*Court Cases (APA)","Judicial Deadline","*Interest Group Opposition")))

Figure1 <- ggplot(Baseline.AME.all, aes(x=AME, y=Var)) + 
  geom_vline(aes(xintercept = 0), size = .25, linetype = "dashed") + 
  geom_pointrange(aes(xmax = upper, xmin = lower), size = .2, position = position_nudge(y = 0), color="black") +
  geom_point(shape=19, size = 2, position = position_nudge(y = 0), color='black') +
  scale_y_discrete(limits=rev) +
  theme(axis.text.y = element_text(hjust=0)) +
  ylab("") +
  xlab("Marginal Effects on Pr(NPRM Rulemaking)") +
  theme_bw(base_size = 14) +
  theme(axis.text= element_text(size=14, colour = "black")) +
  theme(text = element_text(size=14, colour = "black"))
ggsave("Figure1.tiff",width=8.6,height=7,dpi=1200)



# Figure 2
Int.Professionalism.glm <- glm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                                 Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                                 Internal.Joint + Emergency + Contract + 
                                 Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                                 Unalign.President + FirstYear + LastYear +  
                                 Congressional.Opposion + Statutory.Deadline + 
                                 CourtCaseAPA.log + Judicial.Deadline +
                                 Ideo.Mismatch.norm.RCL +
                                 Rule.Novelty:HQ.Pro.per + Routine:HQ.Pro.per + Word.Count.log:HQ.Pro.per + 
                                 Stringency4.log:HQ.Pro.per + Econ.Sig:HQ.Pro.per + SL.Gov:HQ.Pro.per +
                                 External.Joint:HQ.Pro.per +
                                 Agency_Bureau+President, 
                               data=Rule.RCL.Ideology, binomial(link="logit"))
Int.Professionalism.glm.cl <- cluster.vcov(Int.Professionalism.glm, ~Agency)
Int.Professionalism.glm.coef <- broom::tidy(coeftest(Int.Professionalism.glm, Int.Professionalism.glm.cl))

summary(Rule.RCL.Ideology$HQ.Pro.per)
HQ.Pro.per.range=c(HQ.Pro.per=seq(0,82,length.out=83))

Int.AME <- margins(Int.Professionalism.glm, 
                   vcov = cluster.vcov(Int.Professionalism.glm, ~Agency),
                   at = list(HQ.Pro.per = HQ.Pro.per.range), 
                   variables = "Stringency4.log",
                   change = "sd") %>%
  summary()

Figure2 <- ggplot(Int.AME, aes(x = HQ.Pro.per)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .2) +
  geom_line(aes(y = AME)) +
  geom_rug(aes(x = HQ.Pro.per), data = Rule.RCL.Ideology) +
  xlab("Proportion of Professional Bureaucrats (%)") +
  scale_y_continuous(name="Marginal Effects of Rule Stringency on Pr(NPRM Rulemaking)") +
  theme_bw() + 
  theme(axis.text.x = element_text(size = 12, color = "black")) +
  theme(axis.text.y = element_text(size = 12, color = "black")) +
  theme(axis.title.x = element_text(size = 14, color = "black")) +
  theme(axis.title.y = element_text(size = 12, color = "black"))
ggsave("Figure2.tiff",width=7.6,height=5,dpi=1200)


# Figure in Online Appendix C7
levels(Rule.RCL.Ideology$Agency_Bureau)
str(Rule.RCL.Ideology)
Prof.Mean <- Rule.RCL.Ideology %>% 
  dplyr::select(Agency,Bureau,FR.Year,HQ.Pro.per) %>%
  group_by(Agency,Bureau,FR.Year) %>%
  dplyr::summarise(Pro.Mean=mean(HQ.Pro.per,na.rm=T),N=n()) %>%
  group_by(Agency,Bureau) %>%
  dplyr::summarise(Pro.Mean=mean(Pro.Mean,na.rm=T))

pred.pro.SSA = predicts(Int.Professionalism.glm, 
                        "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;5.3;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(146);mode",
                        sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                        sim.count = 10000, set.seed = 123)
pred.pro.FSA = predicts(Int.Professionalism.glm, 
                        "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;8.0;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(9);mode",
                        sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                        sim.count = 10000, set.seed = 123)
pred.pro.AMS = predicts(Int.Professionalism.glm, 
                        "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;11.6;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(6);mode",
                        sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                        sim.count = 10000, set.seed = 123)
pred.pro.FAA = predicts(Int.Professionalism.glm, 
                        "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;25.6;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(104);mode",
                        sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                        sim.count = 10000, set.seed = 123)
pred.pro.NOAA = predicts(Int.Professionalism.glm, 
                         "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;44.7;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(28);mode",
                         sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                         sim.count = 10000, set.seed = 123)
pred.pro.PTO = predicts(Int.Professionalism.glm, 
                        "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;77.5;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(31);mode",
                        sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                        sim.count = 10000, set.seed = 123)
pred.pro.EPA = predicts(Int.Professionalism.glm, 
                        "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;36.9;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(116);mode",
                        sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                        sim.count = 10000, set.seed = 123)
pred.pro.FWS = predicts(Int.Professionalism.glm, 
                        "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;37.0;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(89);mode",
                        sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                        sim.count = 10000, set.seed = 123)
pred.pro.FMCSA = predicts(Int.Professionalism.glm, 
                        "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;18.7;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(106);mode",
                        sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                        sim.count = 10000, set.seed = 123)
pred.pro.VA = predicts(Int.Professionalism.glm, 
                          "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;33.2;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(115);mode",
                          sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                          sim.count = 10000, set.seed = 123)
pred.pro.SEC = predicts(Int.Professionalism.glm, 
                       "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;69.8;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(144);mode",
                       sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                       sim.count = 10000, set.seed = 123)
pred.pro.CFTC = predicts(Int.Professionalism.glm, 
                        "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;52.0;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(2);mode",
                        sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                        sim.count = 10000, set.seed = 123)
pred.pro.NASA = predicts(Int.Professionalism.glm, 
                         "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;35.4;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(128);mode",
                         sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                         sim.count = 10000, set.seed = 123)
pred.pro.IRS = predicts(Int.Professionalism.glm, 
                         "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;32.1;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(99);mode",
                         sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                         sim.count = 10000, set.seed = 123)
pred.pro.FDA = predicts(Int.Professionalism.glm, 
                        "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;68.7;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(47);mode",
                        sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                        sim.count = 10000, set.seed = 123)

pred.pro <- rbind(pred.pro.SSA,pred.pro.AMS,pred.pro.FAA,pred.pro.NOAA,pred.pro.PTO,pred.pro.EPA,
                  pred.pro.FSA,pred.pro.FWS,pred.pro.FMCSA,pred.pro.VA,pred.pro.SEC,pred.pro.CFTC,
                  pred.pro.NASA,pred.pro.IRS,pred.pro.FDA) %>%
  mutate(group=ifelse(Agency_Bureau=="Social Security AdministrationSocial Security Administration","SSA",
               ifelse(Agency_Bureau=="Department of AgricultureAgricultural Marketing Service","AMS",
               ifelse(Agency_Bureau=="Department of TransportationFederal Aviation Administration","FAA",
               ifelse(Agency_Bureau=="Department of CommerceNational Oceanic and Atmospheric Administration","NOAA",
               ifelse(Agency_Bureau=="Department of CommercePatent and Trademark Office","PTO",
               ifelse(Agency_Bureau=="Environmental Protection AgencyEnvironmental Protection Agency","EPA",
               ifelse(Agency_Bureau=="Department of AgricultureFarm Service Agency","FSA",
               ifelse(Agency_Bureau=="Department of the InteriorFish and Wildlife Service","FWS",
               ifelse(Agency_Bureau=="Department of TransportationFederal Motor Carrier Safety Administration","FMCSA",
               ifelse(Agency_Bureau=="Department of Veterans AffairsDepartment of Veterans Affairs","VA",
               ifelse(Agency_Bureau=="Securities and Exchange CommissionSecurities and Exchange Commission","SEC",
               ifelse(Agency_Bureau=="Commodity Futures Trading CommissionCommodity Futures Trading Commission","CFTC",
               ifelse(Agency_Bureau=="National Aeronautics and Space AdministrationNational Aeronautics and Space Administration","NASA",
               ifelse(Agency_Bureau=="Department of the TreasuryInternal Revenue Service","IRS",
               ifelse(Agency_Bureau=="Department of Health and Human ServicesFood and Drug Administration","FDA","XXX")))))))))))))))) %>%
  mutate(group=factor(group,
                      levels=c("SSA","FSA","AMS",
                               "FMCSA","FAA",
                               "IRS","VA","NASA","EPA",
                               "FWS","NOAA","CFTC",
                               "FDA","SEC","PTO")))

label.plot <- data.frame(
  label = c("Avg. Professionals: 5.3","Avg. Professionals: 8.0","Avg. Professionals: 11.6", 
            "Avg. Professionals: 18.7","Avg. Professionals: 25.6", 
            "Avg. Professionals: 32.1", "Avg. Professionals: 33.2","Avg. Professionals: 35.4",
            "Avg. Professionals: 36.9", "Avg. Professionals: 37.0",
            "Avg. Professionals: 44.7", "Avg. Professionals: 52.0",
            "Avg. Professionals: 68.7","Avg. Professionals: 69.8","Avg. Professionals: 77.5"),
  group = c("SSA","FSA","AMS",
            "FMCSA","FAA",
            "IRS","VA","NASA","EPA",
            "FWS","NOAA","CFTC",
            "FDA","SEC","PTO"),
  x = c(6,6,6,6,6,6,6,6,6,6,6,6,6,6,6),
  y = c(0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8)) %>%
  mutate(group=factor(group,
                      levels=c("SSA","FSA","AMS",
                               "FMCSA","FAA",
                               "IRS","VA","NASA","EPA",
                               "FWS","NOAA","CFTC",
                               "FDA","SEC","PTO")))

ggplot(pred.pro, aes(Stringency4.log, mean)) + 
  geom_line(aes(color=group)) +
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=group), alpha=0.15) +
  theme(legend.position="none") + 
  xlab("Stringency") + 
  ylab("Predicted Probability of NPRM Rulemaking") +
  facet_wrap(~group, ncol=5) +
  geom_text(data = label.plot, mapping = aes(x = x, y = y, label = label), size=4) +
  theme(axis.text.x = element_text(size = 12, color = "black")) +
  theme(axis.text.y = element_text(size = 12, color = "black")) +
  theme(axis.title.x = element_text(size = 14, color = "black")) +
  theme(axis.title.y = element_text(size = 14, color = "black"))


pred.pro.EPA.Clinton = predicts(Int.Professionalism.glm, 
                             "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;36.9;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(116);F(2)",
                             sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                             sim.count = 10000, set.seed = 123)
pred.pro.EPA.Bush = predicts(Int.Professionalism.glm, 
                        "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;36.9;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(116);F(1)",
                        sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                        sim.count = 10000, set.seed = 123)
pred.pro.EPA.Obama = predicts(Int.Professionalism.glm, 
                             "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;36.9;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(116);F(3)",
                             sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                             sim.count = 10000, set.seed = 123)
pred.pro.EPA.Trump = predicts(Int.Professionalism.glm, 
                              "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;36.9;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(116);F(4)",
                              sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                              sim.count = 10000, set.seed = 123)

pred.EPA <- rbind(pred.pro.EPA.Clinton,pred.pro.EPA.Bush,pred.pro.EPA.Obama,pred.pro.EPA.Trump) %>%
  mutate(group=ifelse(President=="Bush","Bush",
               ifelse(President=="Obama","Obama",
               ifelse(President=="Trump","Trump",
               ifelse(President=="Clinton","Clinton","XXX"))))) %>%
  mutate(group=factor(group,
                      levels=c("Clinton","Bush","Obama","Trump")))

ggplot(pred.EPA, aes(Stringency4.log, mean)) + 
  geom_line(aes(linetype=group, color=group)) +
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=group), alpha=0.15) +
  theme(legend.position="none") + 
  xlab("Stringency") + 
  ylab("Predicted Probability of NPRM Rulemaking") +
  facet_grid(~group) +
  theme(axis.text.x = element_text(size = 14, color = "black")) +
  theme(axis.text.y = element_text(size = 14, color = "black")) +
  theme(axis.title.x = element_text(size = 14, color = "black")) +
  theme(axis.title.y = element_text(size = 14, color = "black"))


pred.pro.FDA.Clinton = predicts(Int.Professionalism.glm, 
                                "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;68.7;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(47);F(2)",
                                sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                                sim.count = 10000, set.seed = 123)
pred.pro.FDA.Bush = predicts(Int.Professionalism.glm, 
                             "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;68.7;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(47);F(1)",
                             sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                             sim.count = 10000, set.seed = 123)
pred.pro.FDA.Obama = predicts(Int.Professionalism.glm, 
                              "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;68.7;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(47);F(3)",
                              sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                              sim.count = 10000, set.seed = 123)
pred.pro.FDA.Trump = predicts(Int.Professionalism.glm, 
                              "F(1);F(1);mean;0-9,1;F(1);F(1);F(1);F(1);F(1);F(1);mean;68.7;mean;F(1);F(1);F(1);mean;F(1);mean;F(1);mean;F(47);F(4)",
                              sigma=vcovCR(Int.Professionalism.glm, cluster=Rule.RCL.Ideology$Agency, type = "CR2"),
                              sim.count = 10000, set.seed = 123)

pred.FDA <- rbind(pred.pro.FDA.Clinton,pred.pro.FDA.Bush,pred.pro.FDA.Obama,pred.pro.FDA.Trump) %>%
  mutate(group=ifelse(President=="Bush","Bush",
               ifelse(President=="Obama","Obama",
               ifelse(President=="Trump","Trump",
               ifelse(President=="Clinton","Clinton","XXX"))))) %>%
  mutate(group=factor(group,
                      levels=c("Clinton","Bush","Obama","Trump")))

ggplot(pred.FDA, aes(Stringency4.log, mean)) + 
  geom_line(aes(linetype=group, color=group)) +
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=group), alpha=0.15) +
  theme(legend.position="none") + 
  xlab("Stringency") + 
  ylab("Predicted Probability of NPRM Rulemaking") +
  facet_grid(~group) +
  theme(axis.text.x = element_text(size = 14, color = "black")) +
  theme(axis.text.y = element_text(size = 14, color = "black")) +
  theme(axis.title.x = element_text(size = 14, color = "black")) +
  theme(axis.title.y = element_text(size = 14, color = "black"))



# Online Appendix C2: Fixed-Effect Model, Bureau fixed effects
Baseline <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                    Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                    Internal.Joint + Emergency + Contract + 
                    Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                    Unalign.President + FirstYear + LastYear +  
                    Congressional.Opposion + Statutory.Deadline + 
                    CourtCaseAPA.log + Judicial.Deadline +
                    Ideo.Mismatch.norm.RCL |
                    Agency^Bureau + President, 
                  data=Rule.RCL.Ideology, family=binomial(link="logit"),
                  vcov = vcov_cluster("Agency"))
summary(Baseline)

Int.Professionalism <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                               Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                               Internal.Joint + Emergency + Contract + 
                               Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                               Unalign.President + FirstYear + LastYear +  
                               Congressional.Opposion + Statutory.Deadline + 
                               CourtCaseAPA.log + Judicial.Deadline +
                               Ideo.Mismatch.norm.RCL +
                               Rule.Novelty:HQ.Pro.per + Routine:HQ.Pro.per + Word.Count.log:HQ.Pro.per + 
                               Stringency4.log:HQ.Pro.per + Econ.Sig:HQ.Pro.per + SL.Gov:HQ.Pro.per +
                               External.Joint:HQ.Pro.per |
                               Agency^Bureau + President, 
                             data=Rule.RCL.Ideology, family=binomial(link="logit"),
                             vcov = vcov_cluster("Agency"))
summary(Int.Professionalism)

# BH Procedure
Int.Pro.stat <- coeftest(Int.Professionalism)
pval <- Int.Pro.stat[,4]
p.bh <- p.adjust(pval, "fdr")
p.bh05 <- ifelse(p.bh<=0.05, 1, 0)
p.bh10 <- ifelse(p.bh<=0.10, 1, 0)


Int.InterstGroup <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                            Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                            Internal.Joint + Emergency + Contract + 
                            Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                            Unalign.President + FirstYear + LastYear +  
                            Congressional.Opposion + Statutory.Deadline + 
                            CourtCaseAPA.log + Judicial.Deadline +
                            Ideo.Mismatch.norm.RCL +
                            Rule.Novelty:Ideo.Mismatch.norm.RCL + Routine:Ideo.Mismatch.norm.RCL + Word.Count.log:Ideo.Mismatch.norm.RCL + 
                            Stringency4.log:Ideo.Mismatch.norm.RCL + Econ.Sig:Ideo.Mismatch.norm.RCL + SL.Gov:Ideo.Mismatch.norm.RCL +
                            External.Joint:Ideo.Mismatch.norm.RCL |
                            Agency^Bureau + President, 
                          data=Rule.RCL.Ideology, family=binomial(link="logit"),
                          vcov = vcov_cluster("Agency"))
summary(Int.InterstGroup)

# Table C2
dict <- c("New Rule","Routine and Frequent Rule","Rule Complexity",
          "Rule Stringency","Economic significance","Other levels of government",
          "Joint rulemaking - external","Joint rulemaking - internal","Emergency Rule",
          "Contract Rule","Size (employee)","Professionalism",
          "Rulemaking burden","Unaligned president","Presidents’ first year",
          "Presidents’ last year","Congressional opposition","Statutory deadline",
          "Court cases (APA)","Judicial deadline","Interest Group opposition",
          "New Rule X Professionalism","Routine and Frequent Rule X Professionalism",
          "Rule Complexity X Professionalism","Rule Stringency X Professionalism",
          "Economic significance X Professionalism","Other levels of government X Professionalism",
          "Joint rulemaking - external X Professionalism",
          "New Rule X Interest Group opposition","Routine and Frequent Rule X Interest Group opposition",
          "Rule Complexity X Interest Group opposition","Rule Stringency X Interest Group opposition",
          "Economic significance X Interest Group opposition","Other levels of government X Interest Group opposition",
          "Joint rulemaking - external X Interest Group opposition")

texreg::htmlreg(list(Baseline,Int.Professionalism,Int.InterstGroup), 
        file = "Results/Appendix_C2.doc", 
        custom.model.names = c("(1)","(2)","(3)"),
        custom.coef.names = dict,
        digits = 3,
        stars = c(0.001, 0.01, 0.05, 0.1),
        symbol = "&#43",
        inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        custom.gof.rows=list("Bureau Fixed Effects"=c("Yes","Yes","Yes"),
                             "Presidential Administration Fixed Effects"=c("Yes","Yes","Yes")),
        custom.gof.names=c("N",NA,NA,"Deviance","Log Likelihood",NA))



# Online Appendix C3: Fixed-Effect Model, Year-fixed effects 
Baseline.year <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                         Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                    Internal.Joint + Emergency + Contract + 
                    Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                    Unalign.President + 
                    Congressional.Opposion + Statutory.Deadline + 
                    CourtCaseAPA.log + Judicial.Deadline +
                    Ideo.Mismatch.norm.RCL |
                    Agency^Bureau + President + FR.Year, 
                  data=Rule.RCL.Ideology, family=binomial(link="logit"),
                  vcov = vcov_cluster("Agency"))
summary(Baseline.year)

Int.Professionalism.year <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                                    Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                                    Internal.Joint + Emergency + Contract +
                                    Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                                    Unalign.President + 
                                    Congressional.Opposion + Statutory.Deadline + 
                                    CourtCaseAPA.log + Judicial.Deadline +
                                    Ideo.Mismatch.norm.RCL +
                                    Rule.Novelty:HQ.Pro.per + Routine:HQ.Pro.per + Word.Count.log:HQ.Pro.per + 
                                    Stringency4.log:HQ.Pro.per + Econ.Sig:HQ.Pro.per + SL.Gov:HQ.Pro.per +
                                    External.Joint:HQ.Pro.per |
                               Agency^Bureau + President + FR.Year, 
                             data=Rule.RCL.Ideology, family=binomial(link="logit"),
                             vcov = vcov_cluster("Agency"))
summary(Int.Professionalism.year)

Int.InterstGroup.year <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                                 Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                                 Internal.Joint + Emergency + Contract + 
                                 Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                                 Unalign.President + 
                                 Congressional.Opposion + Statutory.Deadline + 
                                 CourtCaseAPA.log + Judicial.Deadline +
                                 Ideo.Mismatch.norm.RCL +
                                 Rule.Novelty:Ideo.Mismatch.norm.RCL + Routine:Ideo.Mismatch.norm.RCL + Word.Count.log:Ideo.Mismatch.norm.RCL + 
                                 Stringency4.log:Ideo.Mismatch.norm.RCL + Econ.Sig:Ideo.Mismatch.norm.RCL + SL.Gov:Ideo.Mismatch.norm.RCL +
                                 External.Joint:Ideo.Mismatch.norm.RCL |
                            Agency^Bureau + President + FR.Year, 
                          data=Rule.RCL.Ideology, family=binomial(link="logit"),
                          vcov = vcov_cluster("Agency"))
summary(Int.InterstGroup.year)

dict.year <- c("New Rule","Routine and Frequent Rule","Rule Complexity",
          "Rule Stringency","Economic significance","Other levels of government",
          "Joint rulemaking - external","Joint rulemaking - internal","Emergency Rule",
          "Contract Rule","Size (employee)","Professionalism",
          "Rulemaking burden","Unaligned president","Congressional opposition","Statutory deadline",
          "Court cases (APA)","Judicial deadline","Interest Group opposition",
          "New Rule X Professionalism","Routine and Frequent Rule X Professionalism",
          "Rule Complexity X Professionalism","Rule Stringency X Professionalism",
          "Economic significance X Professionalism","Other levels of government X Professionalism",
          "Joint rulemaking - external X Professionalism",
          "New Rule X Interest Group opposition","Routine and Frequent Rule X Interest Group opposition",
          "Rule Complexity X Interest Group opposition","Rule Stringency X Interest Group opposition",
          "Economic significance X Interest Group opposition","Other levels of government X Interest Group opposition",
          "Joint rulemaking - external X Interest Group opposition")

texreg::htmlreg(list(Baseline.year,Int.Professionalism.year,Int.InterstGroup.year), 
                file = "Results/Appendix_C3.doc", 
                custom.model.names = c("(1)","(2)","(3)"),
                custom.coef.names = dict.year,
                digits = 3,
                stars = c(0.001, 0.01, 0.05, 0.1),
                symbol = "&#43",
                inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
                head.tag = TRUE, body.tag = TRUE,
                custom.gof.rows=list("Bureau Fixed Effects"=c("Yes","Yes","Yes"),
                                     "Presidential Administration Fixed Effects"=c("Yes","Yes","Yes"),
                                     "Year Fixed Effects"=c("Yes","Yes","Yes")),
                custom.gof.names=c("N",NA,NA,NA,NA,NA,NA))



# Online Appendix C4: Clinton-Lewis ideology estimates (2008)
Baseline.CL <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                       Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                       Internal.Joint + Emergency + Contract + 
                       Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                       Unalign.President.CL + FirstYear + LastYear +  
                       Congressional.Opposion.CL + Statutory.Deadline + 
                       CourtCaseAPA.log + Judicial.Deadline +
                       Ideo.Mismatch.norm.CL |
                       Agency^Bureau + President, 
                       data=Rule.CL.Ideology, family=binomial(link="logit"),
                       vcov = vcov_cluster("Agency"))
summary(Baseline.CL)

Int.Professionalism.CL <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                                  Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                                  Internal.Joint + Emergency + Contract + 
                                  Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                                  Unalign.President.CL + FirstYear + LastYear +  
                                  Congressional.Opposion.CL + Statutory.Deadline + 
                                  CourtCaseAPA.log + Judicial.Deadline +
                                  Ideo.Mismatch.norm.CL +
                                  Rule.Novelty:HQ.Pro.per + Routine:HQ.Pro.per + Word.Count.log:HQ.Pro.per + 
                                  Stringency4.log:HQ.Pro.per + Econ.Sig:HQ.Pro.per + SL.Gov:HQ.Pro.per +
                                  External.Joint:HQ.Pro.per |
                                  Agency^Bureau + President, 
                                  data=Rule.CL.Ideology, family=binomial(link="logit"),
                                  vcov = vcov_cluster("Agency"))
summary(Int.Professionalism.CL)

Int.InterstGroup.CL <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                               Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                               Internal.Joint + Emergency + Contract + 
                               Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                               Unalign.President.CL + FirstYear + LastYear +  
                               Congressional.Opposion.CL + Statutory.Deadline + 
                               CourtCaseAPA.log + Judicial.Deadline +
                               Ideo.Mismatch.norm.CL +
                               Rule.Novelty:Ideo.Mismatch.norm.CL + Routine:Ideo.Mismatch.norm.CL + Word.Count.log:Ideo.Mismatch.norm.CL + 
                               Stringency4.log:Ideo.Mismatch.norm.CL + Econ.Sig:Ideo.Mismatch.norm.CL + SL.Gov:Ideo.Mismatch.norm.CL +
                               External.Joint:Ideo.Mismatch.norm.CL |
                               Agency^Bureau + President, 
                               data=Rule.CL.Ideology, family=binomial(link="logit"),
                               vcov = vcov_cluster("Agency"))
summary(Int.InterstGroup.CL)

texreg::htmlreg(list(Baseline.CL,Int.Professionalism.CL,Int.InterstGroup.CL), 
                file = "Results/Appendix_C4.doc", 
                custom.model.names = c("(1)","(2)","(3)"),
                custom.coef.names = dict,
                digits = 3,
                stars = c(0.001, 0.01, 0.05, 0.1),
                symbol = "&#43",
                inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
                head.tag = TRUE, body.tag = TRUE,
                custom.gof.rows=list("Bureau Fixed Effects"=c("Yes","Yes","Yes"),
                                     "Presidential Administration Fixed Effects"=c("Yes","Yes","Yes")),
                custom.gof.names=c("N",NA,NA,NA,NA,NA))



# Online Appendix C5: 3-digit NAISC Industry Relevance 
Baseline.3digit <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                           Stringency3.log + Econ.Sig + SL.Gov + External.Joint + 
                           Internal.Joint + Emergency + Contract + 
                           Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                           Unalign.President + FirstYear + LastYear +  
                           Congressional.Opposion + Statutory.Deadline + 
                           CourtCaseAPA.log + Judicial.Deadline +
                           Ideo.Mismatch.norm.RCL |
                           Agency^Bureau + President, 
                         data=Rule.RCL.Ideology, family=binomial(link="logit"),
                         vcov = vcov_cluster("Agency"))
summary(Baseline.3digit)

Int.Professionalism.3digit <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                                      Stringency3.log + Econ.Sig + SL.Gov + External.Joint + 
                                      Internal.Joint + Emergency + Contract + 
                                      Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                                      Unalign.President + FirstYear + LastYear +  
                                      Congressional.Opposion + Statutory.Deadline + 
                                      CourtCaseAPA.log + Judicial.Deadline +
                                      Ideo.Mismatch.norm.RCL +
                                      Rule.Novelty:HQ.Pro.per + Routine:HQ.Pro.per + Word.Count.log:HQ.Pro.per + 
                                      Stringency3.log:HQ.Pro.per + Econ.Sig:HQ.Pro.per + SL.Gov:HQ.Pro.per +
                                      External.Joint:HQ.Pro.per |
                                      Agency^Bureau + President, 
                                    data=Rule.RCL.Ideology, family=binomial(link="logit"),
                                    vcov = vcov_cluster("Agency"))
summary(Int.Professionalism.3digit)

Int.InterstGroup.3digit <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                                   Stringency3.log + Econ.Sig + SL.Gov + External.Joint + 
                                   Internal.Joint + Emergency + Contract + 
                                   Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                                   Unalign.President + FirstYear + LastYear +  
                                   Congressional.Opposion + Statutory.Deadline + 
                                   CourtCaseAPA.log + Judicial.Deadline +
                                   Ideo.Mismatch.norm.RCL +
                                   Rule.Novelty:Ideo.Mismatch.norm.RCL + Routine:Ideo.Mismatch.norm.RCL + Word.Count.log:Ideo.Mismatch.norm.RCL + 
                                   Stringency3.log:Ideo.Mismatch.norm.RCL + Econ.Sig:Ideo.Mismatch.norm.RCL + SL.Gov:Ideo.Mismatch.norm.RCL +
                                   External.Joint:Ideo.Mismatch.norm.RCL |
                                   Agency^Bureau + President, 
                                 data=Rule.RCL.Ideology, family=binomial(link="logit"),
                                 vcov = vcov_cluster("Agency"))
summary(Int.InterstGroup.3digit)

texreg::htmlreg(list(Baseline.3digit,Int.Professionalism.3digit,Int.InterstGroup.3digit), 
                file = "Results/Appendix_C5.doc", 
                custom.model.names = c("(1)","(2)","(3)"),
                custom.coef.names = dict,
                digits = 3,
                stars = c(0.001, 0.01, 0.05, 0.1),
                symbol = "&#43",
                inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
                head.tag = TRUE, body.tag = TRUE,
                custom.gof.rows=list("Bureau Fixed Effects"=c("Yes","Yes","Yes"),
                                     "Presidential Administration Fixed Effects"=c("Yes","Yes","Yes")),
                custom.gof.names=c("N",NA,NA,NA,NA,NA))



# Online Appendix C6: PCA + Regression using factor scores
# PCA
Rule.Final.RCL.Factor <- Rule.RCL.Ideology %>%
  dplyr::select(Rule.Novelty,Routine,Word.Count.log,NumberOfPages.log, 
                Stringency4.log,Econ.Sig,SL.Gov,External.Joint) %>%
  mutate_if(is.factor, as.numeric)
  
Uncertainty.PCA <- psych::principal(Rule.Final.RCL.Factor[,c(1:8)],
                                    nfactors=7, rotate="varimax")
Uncertainty.PCA$loadings
Uncertainty.PCA.Score <- data.frame(Uncertainty.PCA$scores)


# Regressions
Rule.Final.RCL.Dim <- cbind(Rule.RCL.Ideology,Uncertainty.PCA.Score) %>%
  rename(Stringency.factor=RC1,Complexity.factor=RC7,
         Econ.factor=RC6,SLgov.factor=RC5,External.factor=RC2,
         Novelty.factor=RC4,Routine.factor=RC3)
colSums(is.na(Rule.Final.RCL.Dim))


Baseline.factor <- feglm(NPRM ~ Novelty.factor + Routine.factor + Complexity.factor + 
                           Stringency.factor + Econ.factor + SLgov.factor + External.factor + 
                           Internal.Joint + Emergency + Contract + 
                           Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                           Unalign.President + FirstYear + LastYear +  
                           Congressional.Opposion + Statutory.Deadline + 
                           CourtCaseAPA.log + Judicial.Deadline +
                           Ideo.Mismatch.norm.RCL |
                           Agency^Bureau + President, 
                         data=Rule.Final.RCL.Dim, family=binomial(link="logit"),
                         vcov = vcov_cluster("Agency"))
summary(Baseline.factor)

Int.Professionalism.factor <- feglm(NPRM ~ Novelty.factor + Routine.factor + Complexity.factor + 
                                      Stringency.factor + Econ.factor + SLgov.factor + External.factor + 
                                      Internal.Joint + Emergency + Contract + 
                                      Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                                      Unalign.President + FirstYear + LastYear +  
                                      Congressional.Opposion + Statutory.Deadline + 
                                      CourtCaseAPA.log + Judicial.Deadline +
                                      Ideo.Mismatch.norm.RCL +
                                      Novelty.factor:HQ.Pro.per + Routine.factor:HQ.Pro.per + Complexity.factor:HQ.Pro.per + 
                                      Stringency.factor:HQ.Pro.per + Econ.factor:HQ.Pro.per + SLgov.factor:HQ.Pro.per +
                                      External.factor:HQ.Pro.per |
                                      Agency^Bureau + President, 
                                    data=Rule.Final.RCL.Dim, family=binomial(link="logit"),
                                    vcov = vcov_cluster("Agency"))
summary(Int.Professionalism.factor)

Int.InterstGroup.factor <- feglm(NPRM ~ Novelty.factor + Routine.factor + Complexity.factor + 
                                   Stringency.factor + Econ.factor + SLgov.factor + External.factor + 
                                   Internal.Joint + Emergency + Contract + 
                                   Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                                   Unalign.President + FirstYear + LastYear +  
                                   Congressional.Opposion + Statutory.Deadline + 
                                   CourtCaseAPA.log + Judicial.Deadline +
                                   Ideo.Mismatch.norm.RCL +
                                   Novelty.factor:Ideo.Mismatch.norm.RCL + Routine.factor:Ideo.Mismatch.norm.RCL + Complexity.factor:Ideo.Mismatch.norm.RCL + 
                                   Stringency.factor:Ideo.Mismatch.norm.RCL + Econ.factor:Ideo.Mismatch.norm.RCL + SLgov.factor:Ideo.Mismatch.norm.RCL +
                                   External.factor:Ideo.Mismatch.norm.RCL |
                                   Agency^Bureau + President, 
                                 data=Rule.Final.RCL.Dim, family=binomial(link="logit"),
                                 vcov = vcov_cluster("Agency"))
summary(Int.InterstGroup.factor)

dict.factor <- c("New Rule Factor","Routine and Frequent Rule Factor","Rule Complexity Factor",
          "Rule Stringency Factor","Economic significance Factor","Other levels of government Factor",
          "Joint rulemaking - external Factor","Joint rulemaking - internal","Emergency Rule",
          "Contract Rule","Size (employee)","Professionalism",
          "Rulemaking burden","Unaligned president","Presidents’ first year",
          "Presidents’ last year","Congressional opposition","Statutory deadline",
          "Court cases (APA)","Judicial deadline","Interest Group opposition",
          "New Rule X Professionalism","Routine and Frequent Rule X Professionalism",
          "Rule Complexity X Professionalism","Rule Stringency X Professionalism",
          "Economic significance X Professionalism","Other levels of government X Professionalism",
          "Joint rulemaking - external X Professionalism",
          "New Rule X Interest Group opposition","Routine and Frequent Rule X Interest Group opposition",
          "Rule Complexity X Interest Group opposition","Rule Stringency X Interest Group opposition",
          "Economic significance X Interest Group opposition","Other levels of government X Interest Group opposition",
          "Joint rulemaking - external X Interest Group opposition")

texreg::htmlreg(list(Baseline.factor,Int.Professionalism.factor,Int.InterstGroup.factor), 
                file = "Results/Appendix_C6.doc", 
                custom.model.names = c("(1)","(2)","(3)"),
                custom.coef.names = dict.factor,
                digits = 3,
                stars = c(0.001, 0.01, 0.05, 0.1),
                symbol = "&#43",
                inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
                head.tag = TRUE, body.tag = TRUE,
                custom.gof.rows=list("Bureau Fixed Effects"=c("Yes","Yes","Yes"),
                                     "Presidential Administration Fixed Effects"=c("Yes","Yes","Yes")),
                custom.gof.names=c("N",NA,NA,NA,NA,NA))



# Online Appendix C8: Presidential administration subsample
Rule.Final.RCL.Clinton <- Rule.RCL.Ideology %>% subset(President=="Clinton")
Rule.Final.RCL.Bush <- Rule.RCL.Ideology %>% subset(President=="Bush")
Rule.Final.RCL.Obama <- Rule.RCL.Ideology %>% subset(President=="Obama")
Rule.Final.RCL.Trump <- Rule.RCL.Ideology %>% subset(President=="Trump")

Baseline.Clinton <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                            Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                            Internal.Joint + Emergency + Contract + 
                            Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                            Congressional.Opposion + Statutory.Deadline + 
                            CourtCaseAPA.log + Judicial.Deadline +
                            Ideo.Mismatch.norm.RCL |
                            Agency^Bureau, 
                       data=Rule.Final.RCL.Clinton, family=binomial(link="logit"),
                       vcov = vcov_cluster("Agency"))
summary(Baseline.Clinton)

Baseline.Bush <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                         Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                         Internal.Joint + Emergency + Contract + 
                         Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                         FirstYear + LastYear +  
                         Congressional.Opposion + Statutory.Deadline + 
                         CourtCaseAPA.log + Judicial.Deadline +
                         Ideo.Mismatch.norm.RCL |
                         Agency^Bureau, 
                  data=Rule.Final.RCL.Bush, family=binomial(link="logit"),
                  vcov = vcov_cluster("Agency"))
summary(Baseline.Bush)

Baseline.Obama <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                          Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                          Internal.Joint + Emergency + Contract + 
                          Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                          FirstYear + LastYear +  
                          Congressional.Opposion + Statutory.Deadline + 
                          CourtCaseAPA.log + Judicial.Deadline +
                          Ideo.Mismatch.norm.RCL |
                          Agency^Bureau, 
                       data=Rule.Final.RCL.Obama, family=binomial(link="logit"),
                       vcov = vcov_cluster("Agency"))
summary(Baseline.Obama)

Baseline.Trump <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                          Stringency4.log + Econ.Sig + SL.Gov + External.Joint + 
                          Internal.Joint + Emergency + Contract + 
                          Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                          FirstYear + LastYear +  
                          Congressional.Opposion + Statutory.Deadline + 
                          CourtCaseAPA.log + Judicial.Deadline +
                          Ideo.Mismatch.norm.RCL |
                          Agency^Bureau, 
                        data=Rule.Final.RCL.Trump, family=binomial(link="logit"),
                        vcov = vcov_cluster("Agency"))
summary(Baseline.Trump)

dict.pres <- c("New Rule","Routine and Frequent Rule","Rule Complexity",
          "Rule Stringency","Economic significance","Other levels of government",
          "Joint rulemaking - external","Joint rulemaking - internal","Emergency Rule",
          "Contract Rule","Size (employee)","Professionalism",
          "Rulemaking burden","Presidents’ first year",
          "Presidents’ last year","Congressional opposition","Statutory deadline",
          "Court cases (APA)","Judicial deadline","Interest Group opposition")

texreg::htmlreg(list(Baseline.Clinton,Baseline.Bush,Baseline.Obama,Baseline.Trump), 
                file = "Results/Appendix_C8.doc", 
                custom.model.names = c("Clinton","Bush","Obama","Trump"),
                custom.coef.names = dict.pres,
                digits = 3,
                stars = c(0.001, 0.01, 0.05, 0.1),
                symbol = "&#43",
                inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
                head.tag = TRUE, body.tag = TRUE,
                custom.gof.rows=list("Bureau Fixed Effects"=c("Yes","Yes","Yes","Yes")),
                custom.gof.names=c("N",NA,NA,NA,NA))



# Online Appendix C9: Separate Measures (Rule Restrictiveness & Scope(Industry relevance))
Baseline.sep <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                    Restrictiveness.log + Industry.Relevance4.log + Econ.Sig + SL.Gov + External.Joint + 
                    Internal.Joint + Emergency + Contract + 
                    Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                    Unalign.President + FirstYear + LastYear +  
                    Congressional.Opposion + Statutory.Deadline + 
                    CourtCaseAPA.log + Judicial.Deadline +
                    Ideo.Mismatch.norm.RCL |
                    Agency^Bureau + President, 
                  data=Rule.RCL.Ideology, family=binomial(link="logit"),
                  vcov = vcov_cluster("Agency"))
summary(Baseline.sep)

Int.Professionalism.sep <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                               Restrictiveness.log + Industry.Relevance4.log + Econ.Sig + SL.Gov + External.Joint + 
                               Internal.Joint + Emergency + Contract + 
                               Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                               Unalign.President + FirstYear + LastYear +  
                               Congressional.Opposion + Statutory.Deadline + 
                               CourtCaseAPA.log + Judicial.Deadline +
                               Ideo.Mismatch.norm.RCL +
                               Rule.Novelty:HQ.Pro.per + Routine:HQ.Pro.per + Word.Count.log:HQ.Pro.per + 
                               Restrictiveness.log:HQ.Pro.per + Industry.Relevance4.log:HQ.Pro.per + Econ.Sig:HQ.Pro.per + SL.Gov:HQ.Pro.per +
                               External.Joint:HQ.Pro.per |
                               Agency^Bureau + President, 
                             data=Rule.RCL.Ideology, family=binomial(link="logit"),
                             vcov = vcov_cluster("Agency"))
summary(Int.Professionalism.sep)

Int.InterstGroup.sep <- feglm(NPRM ~ Rule.Novelty + Routine + Word.Count.log + 
                            Restrictiveness.log + Industry.Relevance4.log + Econ.Sig + SL.Gov + External.Joint + 
                            Internal.Joint + Emergency + Contract + 
                            Tot.Emp.log + HQ.Pro.per + Rule.Burden.log +
                            Unalign.President + FirstYear + LastYear +  
                            Congressional.Opposion + Statutory.Deadline + 
                            CourtCaseAPA.log + Judicial.Deadline +
                            Ideo.Mismatch.norm.RCL +
                            Rule.Novelty:Ideo.Mismatch.norm.RCL + Routine:Ideo.Mismatch.norm.RCL + Word.Count.log:Ideo.Mismatch.norm.RCL + 
                            Restrictiveness.log:Ideo.Mismatch.norm.RCL + Industry.Relevance4.log:Ideo.Mismatch.norm.RCL + Econ.Sig:Ideo.Mismatch.norm.RCL + 
                            SL.Gov:Ideo.Mismatch.norm.RCL + External.Joint:Ideo.Mismatch.norm.RCL |
                            Agency^Bureau + President, 
                          data=Rule.RCL.Ideology, family=binomial(link="logit"),
                          vcov = vcov_cluster("Agency"))
summary(Int.InterstGroup.sep)


dict.sep <- c("New Rule","Routine and Frequent Rule","Rule Complexity",
          "Rule Restrictiveness","Rule Scope","Economic significance","Other levels of government",
          "Joint rulemaking - external","Joint rulemaking - internal","Emergency Rule",
          "Contract Rule","Size (employee)","Professionalism",
          "Rulemaking burden","Unaligned president","Presidents’ first year",
          "Presidents’ last year","Congressional opposition","Statutory deadline",
          "Court cases (APA)","Judicial deadline","Interest Group opposition",
          "New Rule X Professionalism","Routine and Frequent Rule X Professionalism",
          "Rule Complexity X Professionalism","Rule Restrictiveness X Professionalism","Rule Scope X Professionalism",
          "Economic significance X Professionalism","Other levels of government X Professionalism",
          "Joint rulemaking - external X Professionalism",
          "New Rule X Interest Group opposition","Routine and Frequent Rule X Interest Group opposition",
          "Rule Complexity X Interest Group opposition","Rule Restrictiveness X Interest Group opposition","Rule Scope X Interest Group opposition",
          "Economic significance X Interest Group opposition","Other levels of government X Interest Group opposition",
          "Joint rulemaking - external X Interest Group opposition")

texreg::htmlreg(list(Baseline.sep,Int.Professionalism.sep,Int.InterstGroup.sep), 
                file = "Results/Appendix_C9.doc", 
                custom.model.names = c("(1)","(2)","(3)"),
                custom.coef.names = dict.sep,
                digits = 3,
                stars = c(0.001, 0.01, 0.05, 0.1),
                symbol = "&#43",
                inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
                head.tag = TRUE, body.tag = TRUE,
                custom.gof.rows=list("Bureau Fixed Effects"=c("Yes","Yes","Yes"),
                                     "Presidential Administration Fixed Effects"=c("Yes","Yes","Yes")),
                custom.gof.names=c("N",NA,NA,NA,NA,NA))



# Online Appendix C1: Descriptive Statistics
# Rule-level
Descriptive <- Rule.RCL.Ideology %>%
  mutate(Rule.Novelty=as.numeric(as.character(Rule.Novelty)),
         External.Joint=as.numeric(as.character(External.Joint)),
         Internal.Joint=as.numeric(as.character(Internal.Joint)),
         Contract=as.numeric(as.character(Contract)),
         Emergency=as.numeric(as.character(Emergency)),
         Econ.Sig=as.numeric(as.character(Econ.Sig)),
         Routine=as.numeric(as.character(Routine)),
         SL.Gov=as.numeric(as.character(SL.Gov)),
         Statutory.Deadline=as.numeric(as.character(Statutory.Deadline)),
         Judicial.Deadline=as.numeric(as.character(Judicial.Deadline)),
         FirstYear=as.numeric(as.character(FirstYear)),
         LastYear=as.numeric(as.character(LastYear)),
         Unalign.President=as.numeric(as.character(Unalign.President)))
vtable::st(Descriptive,
           vars = c('NPRM',
                    'Rule.Novelty','Word.Count.log','Stringency4.log','External.Joint',
                    'Econ.Sig','Routine','SL.Gov','Internal.Joint','Contract','Emergency',
                    'FirstYear','LastYear','Statutory.Deadline','Judicial.Deadline'), 
           digits=3, fixed.digits = TRUE, numformat = NA)

# Proportions of Interim, Direct, Temporary Final Rules, 2000-2020
N.year <- Rule.RCL.Ideology %>%
  group_by(FR.Year) %>%
  summarise(n=n())

Interim.year <- Rule.RCL.Ideology %>%
  group_by(Interim,FR.Year) %>%
  summarise(n=n()) %>%
  left_join(N.year,by='FR.Year') %>%
  mutate(Interim.per=(n.x/n.y)*100) %>%
  subset(Interim==1)

label.plot <- data.frame(
  label = c("Clinton","Bush","Obama","Trump"),
  x = c(2000,2005,2013,2019),
  y = c(20,20,20,20))

ggplot(Interim.year, aes(x=FR.Year, y=Interim.per)) +
  geom_line() +
  geom_point() +
  ylab("Interim/Direct/Temporary Final Rule (%)") +
  xlab("Year") +
  geom_vline(aes(xintercept = 2001), size = .25, linetype = "dashed") +
  geom_vline(aes(xintercept = 2009), size = .25, linetype = "dashed") +
  geom_vline(aes(xintercept = 2017), size = .25, linetype = "dashed") +
  theme_bw(base_size = 14) +
  scale_y_continuous(limits=c(20,45)) +
  geom_text(data = label.plot, mapping = aes(x = x, y = y, label = label), size=5)


# Bureau-level
Descriptive.Bureau <- Rule.RCL.Ideology %>% 
  dplyr::select(Agency,Bureau,FR.Year,Tot.Emp.log,HQ.Pro.per,Rule.Burden.log,Ideo.Mismatch.norm.RCL) %>%
  group_by(Agency,Bureau,FR.Year) %>%
  dplyr::summarise(Tot.Emp.log.Mean=mean(Tot.Emp.log,na.rm=T),
                   HQ.Pro.per.Mean=mean(HQ.Pro.per,na.rm=T),
                   Rule.Burden.log.Mean=mean(Rule.Burden.log,na.rm=T),
                   Ideo.Mismatch.norm.RCL.Mean=mean(Ideo.Mismatch.norm.RCL,na.rm=T))
vtable::st(Descriptive.Bureau, 
           vars = c('Tot.Emp.log.Mean','HQ.Pro.per.Mean','Rule.Burden.log.Mean','Ideo.Mismatch.norm.RCL.Mean'), 
           digits = 3, fixed.digits = TRUE, numformat = NA)


# Agency-level
Descriptive.Agency <- Descriptive %>%
  dplyr::select(Agency,FR.Year,Unalign.President,Congressional.Opposion,CourtCaseAPA.log) %>%
  group_by(Agency,FR.Year) %>%
  dplyr::summarise(Unalign.President.Mean=mean(Unalign.President,na.rm=T),
                   Congressional.Opposion.Mean=mean(Congressional.Opposion,na.rm=T),
                   CourtCaseAPA.log.Mean=mean(CourtCaseAPA.log,na.rm=T))
vtable::st(Descriptive.Agency, 
           vars = c('Unalign.President.Mean','Congressional.Opposion.Mean','CourtCaseAPA.log.Mean'), 
           digits = 3, fixed.digits = TRUE, numformat = NA)
